{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 第五章　组合优化与模拟退火方法\n",
    "\n",
    "\n",
    "作者：[王何宇](http://person.zju.edu.cn/wangheyu)\n",
    "\n",
    "[浙江大学数学科学学院](http://www.math.zju.edu.cn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这一章我们将不再继续学习如何抽样，而是开始用学到的技术去解决实际优化问题。这里的优化主要指组合优化，尽管我们即将提到的这些方法，对于一些连续优化中的困难情形，如多目标优化等等，也有一定的作用。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 组合优化问题的描述\n",
    "我们知道有很多离散优化问题，是属于NP（Non-deterministic Polynomial ）的，其中最著名的一个可能就是 TSP（Traveling Salesman Problem）问题，它的一种描述如下：给定$n$ 个城市和每两个城市间的距离（无向完全图，有权边），选择一条线路，使得每个城市只经过一次且仅有一次，并最终回到起点。求能完成这个遍历的总路径长度最短。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于完全图，我们用邻接矩阵（adjacency matrix）$D = [d_{ij}]$ 来表示全部城市，其中 $d_{ij}$ 表示城市 $i$ 和 $j$ 之间的距离。显然，$d_{ij}$满足：\n",
    "+ $d_{ij} \\geq 0, 1 \\leq i, j \\leq n$;\n",
    "+ $d_{ii} = 0, 1 \\leq i \\leq n$;\n",
    "+ $d_{ij} = d_{ji}, 1 \\leq i, j \\leq n$;\n",
    "+ $d_{ij} + d_{jk} \\geq d_{ik}, 1 \\leq i, j, k \\leq n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "而问题的一个解可以表达为一个循环排列\n",
    "$$\n",
    "\\pi = (\\pi_1, \\pi_2, \\cdots, \\pi_n),\n",
    "$$\n",
    "表示路径\n",
    "$$\n",
    "\\pi_1 \\to \\pi_2 \\to \\cdots \\to \\pi_n,\n",
    "$$\n",
    "其中 $\\pi_i \\in \\{1, 2, \\cdots, n\\}, i = 1, 2, \\cdots, n$ 表示各个城市的编号，显然对可行解，有 $i \\neq j$，则 $\\pi_i \\neq \\pi_j$。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "若记 $S$ 为全体可行解，则易知\n",
    "$$\n",
    "\\begin{equation}\n",
    "  |S| = \\frac{(n - 1)!}{2}.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "这是一个随 $n$ 增加指数增长的数值，因此对于较大的 $n$ 不太可能用遍历穷举的方式来寻找最优解。这里目标函数可以表示为\n",
    "$$\n",
    "\\begin{equation}\n",
    "  f(\\pi) = \\sum_{i = 1}^n d_{\\pi_i, \\pi_{i + 1}},\n",
    "\\end{equation}\n",
    "$$\n",
    "并约定 $\\pi_{n + 1} = \\pi_1$。\n",
    "\n",
    "于是整个优化问题按组合优化的约定，可以表达为\n",
    "$$\n",
    "\\begin{equation}\n",
    "  \\min f(\\pi_i), \\pi_i \\in S.\n",
    "\\label{eq::TSP_opt}\n",
    "\\end{equation}\n",
    "$$\n",
    "而此优化问题的解则定义为："
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**全局最优解** 称 $\\pi^* \\in S$ 为 TSP 问题的全局最优解，或简称最优解，若 $\\forall \\pi_i \\in S$，有\n",
    "$$\n",
    "f(\\pi^*) \\leq f(\\pi_i).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "作为一个典型的 NP 问题，对每一个可行解（存在的路径），我们都可以迅速得到总路径长度（多项式时间）。然而，求最短路径本身却至今没有什么太好的确定性算法。在这种情况下，随机搜索算法被提到桌面上讨论。然而现在出现新的困难，最优解是一个非常孤立的偶然事件，同时分布上也\n",
    "没有明显的规律（或者我们无法先验地获知它的分布情况），如果完全在 $S$ 全空间随机投点，去碰一个偶然事件，那么收敛速度甚至收敛性本身都很有问题。针对这种未知分布的情况，之前的 Metropolis-Hasting 抽样给我们一定的启示。但首先，我们先观察一下 $S$ 的结构（总是要尽可能获取背景信息）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**$2$-变换邻域** 对于$S$中的一个解\n",
    "$$\n",
    "\\pi = \\{\\pi_1, \\pi_2, \\cdots, \\pi_{p - 1}, \\pi_p, \\pi_{p + 1},\n",
    "\\cdots, \\pi_{q - 1}, \\pi_q, \\pi_{q + 1}, \\cdots, \\pi_n\\},\n",
    "$$\n",
    "称\n",
    "$$\n",
    "N_2(p, q): \\pi \\to \\pi', p < q\n",
    "$$\n",
    "为 $2$-变换，如果\n",
    "$$\n",
    "\\pi'= \\{\\pi_1, \\pi_2, \\cdots, \\pi_{p - 1}, \\pi_q, \\pi_{q - 1},\n",
    "\\cdots, \\pi_{p + 1}, \\pi_p, \\pi_{q + 1}, \\cdots, \\pi_n\\}.\n",
    "$$\n",
    "这里 $\\forall \\pi_i \\in S$，称 $S_i$ 为 $\\pi_i$ 的 $2$-变换邻域，若 $S_i \\subseteq S$，且 $\\forall \\pi_j \\in S_i$，$\\pi_j$ 可由 $\\pi_i$ 经一次 $2$-变换得到。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "类似的，我们还可以定义 $k$-变换，以及 $k$-变换邻域。但注意到 $\\forall i, j \\in S$，我们都可以用至多 $n - 2$ 次 $2$-变换将二者互换（可遍历）。所以从考虑可行解之间关系角度，$2$-变换在很多场合也够用了。有了邻域，我们可以进一步提出：\n",
    "\n",
    "**局部最优解** 记 $\\pi^* \\in S$ 的邻域为 $S^*$，若 $\\forall \\pi_i \\in S^*$，有\n",
    "$$\n",
    "f(\\pi_i) \\geq f(\\pi^*), \n",
    "$$\n",
    "则称 $\\pi^*$ 为 TSP 问题的局部最优解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在我们给出一个用于求解一般组合优化问题的伪代码：\n",
    "```\n",
    "def local_search(i0):\n",
    "    i = i0\n",
    "    while (i不是局部最优解):\n",
    "        从i的邻域Si中选取一个可行解j\n",
    "        if f(j) < f(i):\n",
    "            i = j\n",
    "    return i\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这个算法很没用。它只能提供局部最优解，而 TSP 问题显然要求的是全局最优解或尽可能接近全局最优解。（TSP 算法有实用的离散算法，详见谈之奕老师主讲课程《组合优化》。）我们这里考虑能否使用随机搜索的手段，使这个算法摆脱局部最优解的困扰，但同时我们又要利用局部最优性。因为在没有其他任何信息的时候，局部最优性显然提供了一个邻域内的有效索方向。这里可以考虑的几个关键点是：\n",
    "\n",
    "1. 我们在邻域内每次都做一个随机选取，即总是从 $S_i$ 中随机抽取一个 $j$，作为下一步尝试的可行解。比如具体在 TSP 问题中，可以通过随机抽取 $p$ 和 $q$ 来进行 $2-$变换邻域内的随机抽取；\n",
    "2. 当存在邻域内的可行解 $j$ 满足 $f(j) > f(i)$ 时，也应该有一定概率用 $j$ 去代替$i$，也就是说应该保留一定概率让解往非局部最优的方向行走，这样才能留下脱离局部最优解的可能性；\n",
    "3. 在引入了上述随机搜索之后，原算法的终止判定不再适用，应该相应修改为若当前解 $i$ 在一个指定的步数 $N$（$N >> 1$）都不再变化时，算法终止。\n",
    "\n",
    "这几个问题，在 Metropolis-Hasting 抽样，一一对应都可以解决：\n",
    "1. 给出了实际的邻域选择；\n",
    "2. 我们这里可以假设路径长度在解空间服从 Boltzmann 分布，而引入温度 $T$ 来代表解的稳定程度；\n",
    "3. 就是 Markov 链达到了平稳分布。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Kirkpatrick 在 1983 年将 Metropolis-Hasting 分布抽样方法引入了组合优化，并通过逐步降低温度来得到稳定的近似全局最优解，取得了不错的效果，这一方法被命名为**模拟退火法**。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**定理**  若组合优化问题 $(S, f)$ 的当前解服从平稳分布，则\n",
    "$$\n",
    "\\begin{equation}\n",
    "    \\lim_{t \\to 0}q_i(t) := q_i^* = \\frac{1}{|S_{\\mathrm{opt}}|}\\chi_{\\mathrm{opt}}.\n",
    "  \\end{equation}\n",
    "$$\n",
    "其中 $S_{\\mathrm{opt}}$ 表示最优解集合（不唯一），而 $\\chi_{\\mathrm{opt}}$ 则表示可行解$i$是最优解的特征函数：\n",
    "$$\n",
    "\\begin{equation}\n",
    "    S \\to \\{0, 1\\}, \\chi_{\\mathrm{opt}}(i) = \\left\\{\\begin{array}{ll}\n",
    "    1, & i \\in S_{\\mathrm{opt}},\\\\\n",
    "    0, &\\mbox{其它}.\n",
    "    \\end{array}\\right.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "这个定理表明，如果每个 $t$ 至都能达到平稳分布，那么算法能收敛到全局最优解。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 冷却进度表\n",
    "\n",
    "在实际计算中，模拟退火法有几个重要参数需要决定：初始温度 $t_0$，终止温度 $t_f$，温度衰减函数 $t_k$，以及在每个温度的变换次数或 Markov 链长度 $L_k$。这一节我们面向\n",
    "具体问题来讨论如何确定这些参数。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**初始温度 $t_0$** 确定 $t_0$ 的基本出发点，从 Markov 链出发，我们希望 $t_0$ 基本上是一个无条件准平衡状态。也希望即初始接受率\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "  \\chi_0 = \\frac{\\mbox{接受变换数}}{\\mbox{提出变换数}} \\approx 1.\n",
    "\\end{equation}\n",
    "$$\n",
    "\n",
    "只有这样，我们才能在算法运行初期充分遍历整个可行解空间，发掘潜在的最优解。在这个前提下，根据Metropolis准则\n",
    "$$\n",
    "e^{-\\frac{\\Delta f}{t_0}} \\approx 1,\n",
    "$$\n",
    "知，$t_0$ 应取得较大。或者可以根据需要的 $\\chi_0$ 反向推导 $t_0$。比如如果知道 $\\Delta f$ 的上界为 $100$，且要求 $\\chi_0 = 0.9$，则可取 $t_0 > 949$。但如果缺乏对 $\\Delta f$ 的必要估计，那么也可以用具体的算法过程来控制，比如，Kirkpatrick 1982 提出，可以先确定一个初始 $t_0$，然后做几次尝试变换，若接受率不足，则将 $t_0$ 翻倍，如此循环直到 $t_0$ 满足要求。除此以外还有大量文献讨论这个问题。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**终止温度 $t_f$**\n",
    "一种暴力的做法是限制 $k$ (降温次数)的次数，比如 $k$ 至多 $50$ 次。这种限制主要来自硬件。而从理想角度出发，终止温度恰好和初始温度相反，我们希望它充分小。或者，接受率 $\\chi_f$ 充分小。为此可以分别设置阈值，当 $t_k$ 或 $\\chi_k$ 达到时，就终止算法。而一种更硬核的思路是在连续一个或若干个 Markov 链中，最优解都没有任何变 Kirkpatrick本人就持这种观点。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Markov链长度 $L_k$** 也即在一个 $t_k$ 温度下我们应该迭代多少次。理想的做法当然是迭代到平稳分布，或者至少利用局部信息监测到平稳分布的发生。但这样的算法往往是指数时\n",
    "间的。于是更常用的做法是设置一个多项式的上界 $\\bar{L}$，一旦达到，就不再在此温度继续尝试。比如对 TSP 问题，一个常用的上界是 $\\bar{L} = n^2$，这里 $n$ 指解空间规模 $|S|$ 。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**温度衰减函数 $t_k$** 一个最常用的衰减函数是\n",
    "$$\n",
    "\\begin{equation}\n",
    "  t_{k + 1} = \\alpha \\cdot t_k, k = 0, 1, \\cdots\n",
    "\\end{equation}\n",
    "$$\n",
    "其中 $\\alpha$ 接近 $1$，比如 $\\alpha = 0.95$。（别问为什么！）\n",
    "\n",
    "而稍微有节操一点的做法，比如\n",
    "$$\n",
    "\\begin{equation}\n",
    "  t_k = \\frac{K - k}{K} \\cdot t_0, k = 1, 2, \\cdots, K.\n",
    "\\end{equation}\n",
    "$$\n",
    "这个策略事实上和前面的有限步终止温度是配套的。\n",
    "\n",
    "这里再介绍一个在 TSP 问题中使用良好的降温函数：\n",
    "\n",
    "$$\n",
    "t_k = \\frac{c}{a + \\ln k}, k = 0, 1, 2, \\cdots, K.\n",
    "$$\n",
    "其中 $a$ 和 $c$ 是常数，$c$ 应当充分大以满足 $t_0$ 充分覆盖要求。这种负对数的取法，参考了熵的变化规律，因此在理论和实践中都有较好的表现。\n",
    "\n",
    "以上种种，都是开放问题，在每一个具体问题中，都需要基于上述思路去进行调整。目前在整个领域中，都没有明显的理论最优值可以提供。因此每年还在源源不断地产生大量的论文。大家如果在 SRTP 或者毕业设计中有需要，不妨也考虑一下，不过指导教师建议大家去找更专业的老师。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## TSP 问题的代码实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
